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Large-eddy simulation of a backward facing step 
flow using a least-squares spectral element method 

By Daniel C. Chan 1 AND Rajat Mittal 2 


We report preliminary results obtained from the large eddy simulation of a backward 
facing step at a Reynolds number of 5100. The numerical platform is based on a 
high order Legendre spectral element spatial discretization and a least squares time 
integration scheme. A non-reflective outflow boundary condition is in place to 
minimize the effect of downstream influence. Smagorinsky model with Van Driest 
near wall damping is used for sub-grid scale modeling. Comparisons of mean velocity 
profiles and wall pressure show good agreement with benchmark data. More studies 
are needed to evaluate the sensitivity of this method on numerical parameters before 
it is applied to complex engineering problems. 


1. Introduction 

Many aerospace and commercial products operate in a dynamic flow environment. 
The structural integrity, performance, and development costs of these products are 
affected by the unsteady flowfields they encounter. In rocket propulsion systems, 
dynamic loads are the cause of many life limiting and failure mechanisms. For 
instance, a number of dynamic load related issues manifested themselves during 
the development of the space shuttle main engine, resulting in hundreds of mil- 
lions of dollars of program development costs in terms of hardware redesign and 
testing. Unsteady flows can also be a very effective sound generating mechanism; 
George (1990) states that the aerodynamically generated noise increases approxi- 
mately as velocity to the 6 th power. Therefore, the aerodynamic noise generated 
by vehicles traveling at high speeds can be very annoying to both passengers and 
communities located in the proximity of major highways and railroads. In some 
European countries where trains can travel in excess of 200 MPH, the responsible 
agency has to erect sound walls along the railroads to minimize the effects of noise 
pollution. This requirement can drastically increase the construction and mainte- 
nance costs of a railway system. For passenger cars, unacceptable noise levels inside 
the compartment can have adverse effects on sales. 

In light of the importance in characterizing the dynamic flow environment in 
both aerospace and commercial applications, Rocketdyne has initiated a multi-year 
effort to develop a general purpose computational fluid dynamics based analysis 
system for dynamic load prediction. This system will provide high-fidelity pre- 
dictive capability through the development of a novel numerical algorithm and 
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utilization of distributed parallel computing. The numerical algorithm is a high 
order spectral method which provides the unique capability to accurately model 
complex geometries and rapidly varying flowfields. Parallel computing provides 
the necessary memory capacity and speed required for large scale computations. 
All these features have been incorporated in the Rocketdyne Unstructured Implicit 
Flow (Uni Flo) solver. The UniFlo code is capable of performing a hierarchy of fluid 
dynamic analyses including direct numerical simulation (DNS), large eddy simula- 
tion (LES) and Reynolds average Navier-Stokes solution (RANS). Only DNS and 
LES can provide time accurate information that is needed for unsteady turbulent 
simulations. LES models flow features that are not directly captured by the grid 
resolution employed. This technique is also known as subgrid scale (SGS) modeling. 
The LES approach (vs. DNS) can relax the requirement on grid resolution that is 
normally very demanding for turbulent flow simulations making it an effective tool 
for engineering analyses. However, one also has to be concerned with the numerical 
errors that increase as the grid is coarsened. If not controlled properly, these errors 
can overwhelm the advantage offered by LES. Therefore, the purpose of this work 
is to first evaluate the numerical accuracy of UniFlo in predicting time dependent 
flows. Once this is accomplished, we then assess the capability of the Smagorinsky 
SGS model in predicting turbulent flow. The backward facing step configuration is 
chosen as the benchmark case since it mimics the flowfield in a rocket engine com- 
bustor and existing numerical and experimental data are available for comparison. 

In what follows, we describe the numerical method, boundary condition and 
SGS model employed by UniFlo. Numerical results demonstrating accuracy of the 
method and effectiveness of the Smagorinsky model are also provided. 


2. Numerical method 

The Navier-Stokes equations are written as a first order system and can be repre- 
sented as Cu — f in a domain Q C $i nrf which is subjected to the condition Bu = g 
along a piecewise smooth boundary T. £ is a first-order partial differential operator: 


n d 


du 




rid = 2 or 3, depending on the spatial dimensions, x\s are the Cartesian coordinates, 
u has a length n, where n is the number of dependent variables, / is the forcing 
function, and both B and g describe the appropriate boundary conditions. A's are 
m x n matrices which describe the characteristics of the system of equations being 
solved. The idea behind the least squares spectral element method (LSSEM) is to 
minimize the residual 

R = Cu — f 

in a least squares sense within the domain of interest and construct the functional 
as 

I(u) = l - II Cu - f\\l = (j Cu - f, Cu - /) 
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By setting 61 = 0 and 6u = u>, one can reduce the problem to 

(. Cw,Cu ) = ^£w,cfj w € 5 

where, 5 = {u € /fo(S7);.Bff = g on T}, and is the Sobolev space with a 
compact support. For incompressible viscous flows, the working variables are ve- 
locity, pressure, and vorticity. By using this system of equations, one can employ 
any of the C° functions to approximate the spatial variation of the dependent vari- 
ables. UniFlo employs isoparametric mapping to transform the governing equations 
from the Cartesian coordinate system to a generalized coordinate system where the 
spatial discretization is performed. The domain of interest is divided into a set 
of non-overlapping elements and within each element, basis function derived from 
Legendre polynomials is used for spatial discretization. The spatial accuracy de- 
pends on the choice of the order of Legendre polynomial basis function and can vary 
from element to element. This approach, also known as spectral element, has been 
formulated by R0nquist and Patera (1987). LSSEM uses a common interpolating 
function to approximate all of the dependent variables. Even with the presence of 
the convective terms, the resulting set of algebraic equations are positive definite 
and symmetric. LSSEM maintains a tight coupling among all of the governing 
equations and provides a set of well-defined boundary conditions that are consis- 
tent with flow physics and mathematical constraints. It does not require any user 
defined artificial damping factor to maintain numerical stability. To maintain high 
spatial accuracy at the domain boundary, UniFlo does not need special treatment 
such as the utilization of ghost points. The convective terms axe linearized with the 
Newton-Raphson procedure so that the spatial derivatives can be discretized im- 
plicitly. Sub-iterations are required at each time step for the purpose of minimizing 
the effect of linearization errors. For most problems, the residual can be reduced 
by four orders of magnitude in less than three iterations. The accuracy is second 
order in time with the application of a backward differencing scheme. For instance, 
the temporal derivative of the velocity component, u, can be discretized as 

du u 9 ~ 1 — 4u 3 + 3u 9 ^ 1 
~dt ~ 2 6 t 

where superscripts represent different time levels. The resulting algebraic equations 
are solved by the conjugate gradient method with Jacobi preconditioning. The 
structure of the coefficient matrix is completely arbitrary and the solution procedure 
does not rely on any pre-defined order. More details of this method is given by 
Chan (1996). 

The boundary conditions are: (1) specified velocity at the inlet, (2) no slip along 
solid walls, (3) stress free and vanishing normal velocity component along the plane 
of symmetry and (4) Tree boundary’ along an outflow plane. For a Cartesian grid, 
stress free condition is imposed by setting the horizontal vorticity components to 
zero. Points located on a Tree boundary’ are treated as unknowns and solved 
directly. 
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FIGURE 1. Streamlines behind a backward facing step at Re = 389; top half: 
result for long domain; bottom half: result for truncated domain. 


For turbulent flows, we relate the subgrid scale stresses to the strain rate of the 
resolved velocity field via Boussinesq approximation. The diffusion term of the 
Navier-Stokes equations then becomes 


1 du k dv, 

{ Re +Ut),}k d Xj d Xj S,i 


where u t is the eddy viscosity, S tJ is the strain rate, and u is the vorticity. The value 
of eijk is equal to zero unless each of the number 1,2, and 3 occurs as a subscript. 
Furthermore, e,j* is equal to 1 if the order of subscripts is cyclic, it becomes -1 if 
the order of subscripts is not cyclic. The eddy viscosity is computed as 


= (C,A) 2 \S, } \f a 


A = ( 6x6y6z )^ 3 

where C s = 0.1 and f s is the Van Driest damping function defined as 

/« = l.O-exp(-l-) 

In reality the value of C 9 is not constant and can change in time and space. Near a 
corner, 8 + is determined with the shortest normal distance from the adjacent walls. 
This procedure is somewhat ad hoc and is problem dependent. 
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FIGURE 2. Predicted profiles at an axial distance of 5 inlet heights behind a 

backward-facing step with Re = 389, results obtained with long domain 

results obtained with truncated domain, (a) axial velocity profile and (b) 

vorticity profile. 

3. Numerical results and discussion 

To demonstrate the effectiveness of the current outflow boundary condition, we 
apply it to compute the laminar flow behind a backward facing step studied ex- 
perimentally by Armaly et al. (1983). The Reynolds number, based on the inlet 
height and average velocity, is 389. The ratio between the inlet and step heights 
is 0.94. Flow separates behind the step and reattaches at an axial distance that is 
equal to about eight step heights from the plane of expansion. Two exit domains, 
one long and one short, are used. For the long domain case, the axial length be- 
hind the step is 17, the flow has room to reattach after separation and recover to a 
fully- developed flow; therefore, the downstream influence on the flowfield near the 
step is small, and for comparison we can use the predicted profiles as the baseline. 
For the short domain case, the outflow plane, which cuts through the separated 
region, is located at 5 inlet heights behind the step. Because of this, accuracy 
of the predicted profiles is strongly influenced by the outflow boundary condition. 
For time dependent turbulent flow, this situation is similar to having an eddy pass 
across an outflow boundary. A parabolic profile is imposed along the inlet plane 
which is located at 2 inlet heights upstream of the expansion. Figure 1 shows the 
grid systems and streamlines predicted by UniFlo for both the short and long do- 
mains. In both cases, 5 collocation points are placed within each element. The 
total number of elements is 72 for the long domain and 36 for the short domain. 
The flow pattern is almost the same in both cases. For the short domain case, 
having a reverse flow on part of the outflow boundary does not present numerical 
convergence problem, and this further demonstrates the robustness of the current 
numerical method and outflow boundary condition. The predicted reattachment is 
8.0 times the inlet height and is in good agreement with the test data. Armaly et 
al. also reports that at Re — 389, the flow begins to separate from the upper wall 
and becomes three-dimensional, but the separation region is so small that its size 
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FIGURE 3. Predicted profiles behind a backward-facing step with Re = 800; 

Gartling’s results, o 5 th order, A 6 th order, and □ 7 th order; (a) axial location 

of 7 and (b) axial location of 15. 

could not be measured. This phenomena is correctly predicted by UniFlo. Figure 2 
shows the axial velocity and vorticity profiles at an axial location of 5 inlet heights 
behind the step. The trend in both cases is identical, with only less than 10 percent 
discrepancy on the magnitude. 

The next test case is due to Gartling (1990) and Gresho et al. (1993). The pur- 
pose of this exercise is to answer some of the questions raised by Gresho et al as 
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FIGURE 4. Predicted wall vorticity distribution for a backward-facing step with 
Re = 800, o 5 th order, A 6 <fc order, and □ 7 th order, (a) upper wall and (b) lower 
wall. 


to whether spectral methods can handle flow geometries with a sharp corner and 
predict the correct flow behavior. Through careful numerical studies and stability 
analysis, they conclude that at a Reynolds number of 800, the flow behind a back- 
ward facing step with 1:2 expansion ratio is indeed steady. With this in mind, we 
first perform the simulation as a steady state problem by turning off the transient 
terms in the Navier-Stokes equations. The rectangular flow domain is 17 units long 
and 1 unit high. The flow enters the domain along the top half portion of the left 
boundary with a parabolic profile. The Reynolds number based on the step height 
and mean velocity is 800. Figure 4 shows the grid skeleton employed; there are 
4 elements in the vertical direction and 11 elements in the streamwise direction. 
Within each element, we apply 5* h , 6 th , and 7 th order polynomials, respectively, in 
each of the two directions. Figure 3 shows the comparison between the predicted 
profiles and benchmark data at two different streamwise locations. All except the 
vertical velocity profile at the axial location of 7 show an excellent agreement with 
the benchmark data of Gartling. Figure 4 shows the vorticity distribution, which 
is proportional to shear stress, along the bottom and top boundaries. By examin- 
ing these plots, one can determine both the separation and reattachment points. 
Along the lower wall, UniFlo predicts a reattachment length of 6.1, whereas along 
the upper wall, it predicts a separation at the streamwise location of 4.8 and a 
reattachment at the streamwise location of 10.5. These predictions are in excellent 
agreement with the benchmark data. These results also indicate that for steady 
flow computation, numerical error incurred from using an under-resolving grid is 
very localized. 

We then compute the same problem by treating it as an unsteady flow. Initially, 
the flow is stagnant inside the domain. Figure 5 shows the temporal evolution of 
the streamlines for the case where 6 th order polynomials are used inside each el- 
ement. Overall this grid resolution produces satisfactory results for steady state 
calculation, however, this is not the case for time accurate simulation. A transient 
process, which involves a sequence of vortex shedding, takes place along the upper 
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FIGURE 5. Streamlines showing the time evolution of the flowfield behind a 
backward facing step at a Reynolds number of 800; computation performed with a 
11x4 grid and 6 th order Legendre polynomials; final state is a temporally periodic 
flow; from top to bottom: time=10, 20, 30, 50, 80, 100, and 140. 

wall at the streamwise location where the steady state result show a discrepancy 
in the vertical velocity profile prediction. This result demonstrates that numeri- 
cal error that develops in a small region can grow over time and contaminate the 
entire flowfield. We then refine the grid by increasing the? number of elements in 
the streamwise direction to 18 while maintaining 6 th order polynomials in each ele- 
ment. The result shown in Fig. 6 indicates that the initial transient flow features 
decay rapidly in time and the flow evolves asymptotically towards a steady state. 
This prediction agrees with the finding of Gresho v.t al. It is apparent that the 
transient flow predicted above is a numerical artifact. Unfortunately, the flow fea- 
tures generated by this numerical error look so real, making them difficult to detect. 
Therefore, for unsteady flow simulation one must perforin grid dependence study 
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FIGURE 6. Streamlines showing the time evolution of the flowfield behind a 
backward facing step at a Reynolds number of 800; computation performed with a 
18 x 4 grid and order Legendre polynomials; asymptotic state is steady; from 
top to bottom: time=10, 20, 30, 50, 80, 100, and 140. 

before attempting to explain the underlying flow physics. 

Having addressed some of the relevant numerical issues, we then use Uni Flo to 
simulate the three-dimensional backward-facing step configuration where experi- 
mental data of Jovic and Driver (1994), DNS data of Le and Moin (1994), and LES 
data of Akselvoll and Moin (1995) are available for comparison. The grid system 
employed is shown in Fig. 7. There are 13 elements in the streamwise direction, 6 
elements in the vertical direction and 6 elements in the spanwise direction. Within 
each element, order Legendre polynomials are used in each of the three direc- 
tions. The expansion ratio is 5:6. The geometry is scaled with the step height, H . 
The inlet plane is located at 5 H upstream of the expansion, and the outflow plane is 
located at 17 H downstream of the expansion plane. The spanwise width is 4 H and 
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FIGURE 7. Grid system for the three-dimensional backward facing step, top: 
through flow plane, bottom: cross-sectional plane. 


a periodic boundary condition is imposed in this direction. Stress free condition is 
imposed along the top boundary, and no slip condition is imposed along the bottom 
wall. At the inlet plane, we take the time dependent turbulent boundary layer pro- 
files computed by Akselvoll and Moin and interpolate them onto the current grid. 
The freestream velocity is taken to be one, and the time it takes for the flow to 
travel one step height is also one. Since our implicit flow solver is not restricted 
by the CFL condition for numerical stability, we can take a larger time step size 
of 0.1, which is five times higher than that employed by Akselvoll and Moin. As a 
result, each through flow takes 220 time steps. Time average quantities are collected 
after the flow has evolved through the domain 5 times. For comparison, we further 
average the data in the spanwise direction and show them in Fig. 8. The predicted 
wall pressure distribution is in good agreement with the experimental data inside 
the recirculating region behind the step. However, in the recovery region, all the 
numerical methods, including UniFlo, predict a faster recovery rate than that of the 
experimental measurements. The velocity profile's are compared to the DNS data 
of Le and Moin. The agreement, is good for all five stream wise locations. 
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FIGURE 8. Time and spanwise average flow quantities predicted by UniFlo for a 
three-dimensional backward facing step at Re — 5100, top: pressure coefficient along 

bottom wall, UniFlo prediction, DNS result of Le &; Moin LES 

result of Akselvoll Sz Moin, and o experimental result of Jovic &z Driver; bottom: 

axial velocity profiles at selected streamwise locations; UniFlo, o DNS result 

of Le $z Moin. 

4. Summary 

We have demonstrated that a spectral based flow solver can be used to simulate 
the flow behind a backward-facing configuration. The weak singularity located 
at the corner does not present numerical problem to the least squares method. 
Numerical error can generate unsteady flow phenomena that could be mistaken as 
‘real’ flow physics. Therefore, grid dependence study is paramount (more so than 
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the steady state flow calculation) in unsteady flow simulation. The preliminary 
results obtained from the LES show good agreement with both experimental data 
and numerical data. Further studies are needed in order to understand the role of 
the subgrid scale model in these simulations. The Smagorinsky model with Van 
Driest wall damping is, however, difficult to implement for complex geometries. 
Future work will include the implementation of dynamic models that do not require 
wall damping function and user specified model constant. 
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